3D Euler equations and ideal MHD mapped to regular systems: 
probing the finite-time blowup hypothesis 
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We prove by an explicit construction that solutions to incompressible 3D Euler equa- 
tions defined in the periodic cube = [0,L]^ can be mapped bijectively to a new 
system of equations whose solutions are globally regular. We establish that the usual 
Beale-Kato-Majda criterion for finite-time singularity (or blowup) of a solution to the 
3D Euler system is equivalent to a condition on the corresponding regular solution of 
the new system. In the hypothetical case of Euler finite-time singularity, we provide an 
explicit formula for the blowup time in terms of the regular solution of the new system. 
O The new system is amenable to being integrated numerically using similar methods as in 

— . Euler equations. We propose a method to simulate numerically the new regular system 

psj and describe how to use this to draw robust and reliable conclusions on the finite-time 

O singularity problem of Euler equations, based on the conservation of quantities directly 






related to energy and circulation. The method of mapping to a regular system can be 

extended to any fluid equation that admits a Beale-Kato-Majda type of theorem, e.g. 3D 

^^ Navier-Stokes, 2D and 3D magnetohydrodynamics, and ID inviscid Burgers. We discuss 

I— I briefly the case of 2D ideal magnetohydrodynamics. In order to illustrate the usefulness 

►j^ of the mapping, we provide a thorough comparison of the analytical solution versus the 

'^ numerical solution in the case of ID inviscid Burgers equation. 

i 

^^ PACS numbers: 47.10.A-, 47.ll.-j, 47.15. ki, 47.65.-d 
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I. INTRODUCTION 

One of the most important unsolved problems in Mathematics entails a simple question: are 
solutions to the three-dimensional Euler equations globally regular or do they blow up in a finite 
time? The analogous question for Navier-Stokes equations, also unsolved, corresponds to one 
of the famous Millennium Prize problems [Ij. The analytical and numerical methods appearing 
in the scientific literature to solve these equations have been highly transferrable to real-life 
problems, such as high-Reynolds number turbulence and vortex reconnection in Navier-Stokes, 
magnetic reconnection in magnetohydrodynamics, and extreme events in the atmosphere, to 
mention a few. 

While several conclusive results are available for Euler and Navier-Stokes in two dimensions, 
three dimensions with axial symmetry and three dimensions with helical symmetry ^-S\, attempts 
to understand the regularity of three-dimensional solutions have only reached local and/or con- 
ditional analytical results [9-28J. Since the advent of powerful computers circa 1980, numerical 
researchers have also contributed with competing yes/no conclusions regarding finite-time sin- 
gularity, mainly in three-dimensional Euler, which is the focus of the present paper [22H55]. 

In a 3D Euler numerical simulation the spatial distribution of vorticity tends to get localised in 
structures that become increasingly sharp with time. This entails a finite-time loss of resolution, 
not necessarily due to a true finite-time singularity of the solutions, but rather to a finite amount 
of memory available for a computer simulation. One can identify two main drawbacks in current 
and previous state-of-the-art numerical attempts on the Euler singularity problem. On the one 
hand, it is not known analytically whether a given initial condition can give rise to a finite-time 
singularity, hence a numerical solution of 3D Euler equations is inherently "blind" to any potential 
finite-time singularity that may be encountered: numerical dissipation and loss of resolution may 
certainly "shield" the singularity. On the other hand, instability and numerical error could give 
rise to a "fake" finite-time singularity, typically associated with a lack of resolution at late stages 
of a simulation. 

This observation motivated the current research: our main analytical result is that there is a 
bijective mapping from Euler equations (along with velocity fields) to a new system of equations, 
hereby called mapped, whose solutions are globally regular. The mapped system is amenable to 
direct numerical integration using the same methods as in Euler equations, with the important 
advantage that the solutions of the mapped system are regular by definition. In this way, and 
for the first time ever, reliable and robust conclusions on the problem of finite-time singularity of 
Euler equations could be drawn, by analysing in post-processing the data from a careful numerical 
integration of the mapped system. 

In Section [M] we define the bijective mapping from 3D Euler to a regular system, provide 
the proofs of global regularity, and show how to use the solution of the regular system to 
probe the hypothesis of finite-time singularity of the 3D Euler equations. In Section [ITI] we 
propose a numerical method to implement these ideas and discuss the advantages and potential 



challenges of the new approach. In Section [lV]we show that other equations of interest in physics 
and mathematics can be mapped to regular systems, provided they admit a Beale-Kato-Majda 
(BKM) type of theorem (e.g., Navier-Stokes, magnetohydrodynamics and inviscid Burgers), and 
we construct an explicit mapping from the 2D ideal magnetohydrodynamics equations to a 
regular system. In Section [V] we illustrate the usefulness of the mapping to a regular system, 
by considering a simple equation admitting a BKM type of theorem: the ID inviscid Burgers 
equation. We compare thoroughly the analytical solution for the supremum norm of velocity 
gradient, with the numerical data from both direct numerical simulation of inviscid Burgers and 
numerical simulation of the mapped system, and conclude that the integration of the regular 
mapped system is far more accurate and converges better than the original system. Finally, we 



present concluding remarks in Section VI 



II. EULER EQUATIONS: PREVIOUS VS. NOVEL APPROACH 

A. Euler equations in a periodic cuboid 

The Euler fluid equations describe the evolution of an incompressible, unit mass density flow 
with velocity field u{x,y,z,t) G M^ defined for {x,y,z) G M? and in a time interval t G [0,T): 

^ + u ■ Vu = -Vp , V ■ u = 0. (1) 

at 

Periodic boundary conditions are assumed for the velocity field and pressure with a basic pe- 
riodicity domain Vt = [0,L]^. That is, for any {x,y,z) G M? we have u(x + L,y,z,t) = 
u(x, y + L,z,t) = u{x, y,z + L,t) = u(x, y, z, t),\/t G [0, T) and similarly for p. The basic 
periodicity domain is by definition the smallest domain with such periodicity property. 

Here, and throughout this proposal, T denotes a generic time so that u G C([0,T);/7*) fl 
C^([0, T); if '*~^) , s > 3 , so in particular the Sobolev norms of the velocity field are bounded: 



|u(-,t)||^. = \}^{l + \k\y\u{k,t)\-'\ <cs, tG[0,T), .>3, (2) 

\kez3 / 

where Cg are some constants and u(k, t) are the Fourier coefficients of the velocity field. 

B. Previous analytical and numerical results: vorticity and the finite-time singularity 

hypothesis 

In order to gain some insight on the problem of finite-time blowup, the most common approach 
is to monitor from numerical simulations a number of quantities pertaining to the theory. The 
main quantity studied since early numerical simulations (ca. 1987) is the vorticity supremum 
norm ||a;(-,t)||oo = supxgf7|<^(x, t)| , where lj = V x u is the vorticity field. According to a 
theorem by Beale, Kato and Majda (BKM) |T3], the assumed regularity of the velocity field can 
be extended up to and including the time T if and only if t{T) = J^ ||a;(-, t)||oo dt < oo. By 
'regular up to and including the time T' we mean uG C{[0,T]; H') n C^([0,T]; H'^^), s > 3. 

The hypothesis of finite-time singularity (blowup) states that there is a 'singularity time' 
T G (0, oo), such that t(T) = cxd. It is an open problem to establish analytically the validity 
of this hypothesis. Consequently, it has been an important aim of many numerical simulations 
of Euler equations ([I]), published since BKM theorem emerged, to conclude yes, no or maybe 
on the validity of this hypothesis. Unfortunately, direct numerical simulations of Euler equations 
entail an inherent uncertainty about the time of a hypothetical singularity. Therefore, current 
numerical attempts on the Euler singularity problem cannot give 100% reliable conclusions in 
favour of or against a finite-time singularity. A fresh approach is urgently needed to continue 
along this research path. 

C. The novel approach 

The main goal of this paper is to provide such a fresh approach to study the validity of the 
hypothesis of finite-time singularity. This is based on a bijective mapping which maps the Euler 
equations to a system whose solutions are globally regular. The mapping from 'old variables' 
(t,u(x, t)) to 'mapped variables' (r, Urnap(x, r)) is defined by: 

T{t)= [ \\u:{-,t')\\^dt', u^ap(x,r)= "^""'Y ■ (3) 

Jo \\^{-,tj\\oo 



The mapped time r(t) is a strictly monotonically increasing function of t and the mapped 
velocity is well defined due to the following result: 

Lemma 1. Assume the initial vorticity is not identically zero on Vl. Then the vorticity supremum 
norm is positive: ||lj(-, t)||oo > 0,Vt G [0,T). D 

Proof. Using the well-known fact that the zeroes of vorticity are preserved by the flow for as 
long as the vorticity is defined, we get: ||uj(-, ti)||oo = for some ti G [0,T) if and only 
if ||li;(-, t)||oo = 0, Vt G [0,T). Since the initial vorticity is not identically zero, the Lemma 
follows. D 

Moreover, the characteristics of each system are mapped accordingly. Let X(Xo,t) and 
XmaplXcr) be the respective solutions of 

— -X(Xo, t) = u(X(Xo, t),t) , — Xmap(Xo, t) = Umap(Xmap(Xo, t), t) , 

with initial conditions Xinap(Xo,0) = X(Xo,0) = Xq . Then, using equations ([S]) we get 
X^ap(Xo,r(t)) = X(Xo,t). 

The Euler equations ([I]) are equivalent to the system: 

—^^ + Umap ■ VUmap = -Vpmap - /3(r) U^ap , V ■ U^ap = , (4) 

where all mapped fields are to be evaluated at (x, r) unless explicitly stated, and the coefficient 
/3(r) (damping if and only if J^ |uinap(x, r)p (Px decreases in time) is defined by 

/3(r) = a;^ap(Y(r), r) • (Vu^ap(Y(r), r)) • a;,,ap(Y(r), r) , (5) 

where LJinap(x, r) = V x Uinap(x, r) is the mapped vorticity and Y(r) is the position of mapped 
vorticity maximum, with |a;map(x, r)| < ||u;map(-, ^)||oo = |'*^map(Y(r), r)| = 1 Vr, Vx G fi. 
The existence of Y(r) and /3(r) is a consequence of the assumed regularity of the velocity 
field: u G C{[0,T);H') n C\[0,T); H'-'^) , s > 3, along with Lemma 1 and a well-known 
Sobolev lemma: 

\\D^F{-)\\^ < c, ||F(-)|U.+P , Vp G Z+ U {0}, s > 3/2 , 

valid for 3D vector fields F defined on fi, where D^ is any combination of spatial derivatives of 
combined order p. 

In order to be able to treat Y(r) as a function, a selection mechanism is required. Apart from 
mirror images -which do not pose any problem for the selection mechanism-, in the generic case 
Y(r) can be taken as a piecewise continuous function, with jumps at times when two competing 
local spatial maxima of |ci;map(x, r)|, located at different places, coincide in value. The full 
treatment of this generic case of discontinuous Y(r) is possible and straightforward. But, for the 
sake of this paper and to fix ideas, we assume from here on that Y(r) is continuous. We remark 
that the main results discussed in this paper (Theorem 2 and Corollary 3) can be extended 
with slight modifications to the general case of discontinuous Y(r). It is worth mentioning 
that several initial conditions considered as candidates for finite-time singularity [321 [341 |4U1 I5l] 
belong to the simplified case of continuous Y(r). 

The new system (|4])-(|5]) is globally regular: 

Theorem 2. The solution of the new system ^-^ is regular for all finite values of the mapped 
time < r < oo. D 



Proof. Let us suppose r G [0, tq] with < tq < oo fixed. We obtain, using equation ([S]), 

f \\Lj{;t")\\^dt"<To<00, Vt'G[0,t(ro)], 

where t^r) is the inverse map from new to original time. According to BKM theorem 
this implies that in original variables all Sobolev norms of velocity ||u(-,t)||i:/ii are bounded for 
all times t E [0,t(ro)]. Also, by virtue of Lemma 1 the vorticity modulus supremum norm is 
positive: ||uj(-, t)||oo > 0,t G [0,t(ro)]. Therefore we can transform to the mapped variables 
and use equation ([2]) to conclude that all Sobolev norms of mapped velocity ||umap(-7 t)||//^' = 
||u(-,t(r))||//s/||Lj(^, t(r))||oo are bounded for r G [0,ro]. Since tq is arbitrary the result 
Umap ^ C{[0,oo); H^) , s > 3 is established. Along similar lines we can show that u^ap £ 
C^{[0,oo);H''~'^). The proof is omitted here. D 

Theorem 2 holds even if the Euler system ([I]) has a finite-time singularity. This is very useful 
from a computational point of view: one is sure that a careful numerical integration of the new 
system should give rise to a regular solution, and this alone leads to a more reliable test of the 
validity of the finite-time singularity hypothesis. This is in contrast with numerical simulations of 
the original Euler equations, where numerical errors and lack of resolution may fail to find true 
finite-time singularities or may give rise to spurious ones, making any conclusions unreliable. 

D. A robust criterion for finite-time singularity of Euler system 

In the original variables there are two types of constants of motion: the total kinetic energy 
of the fluid and the circulations of velocity field along selected closed contours. For simplicity we 
focus on the energy E = ^||u(-, t)||^2 = |/n|u(x, t)p (Px, defined in terms of the L^-norm of 
velocity. In pseudospectral Euler simulations, constancy of energy is very robust and it holds even 
when small-scale errors pervade the system. In a numerical simulation of the mapped system 
(|4])-([5]) we do not have direct access to the Euler energy but we can define the mapped energy: 

1 If 

^map(T) = 2l|Umap(-,T)|||2 = " / |Umap(x,r)P d^X, (6) 

which is not constant. From ([S]) we deduce that the supremum norm of the original vorticity 
||lj(-, t(r))||oo can be obtained in terms of the Euler energy and mapped energy: 



l^(-,i(r))||oo = W;^;;;^' (7) 

where the implicit inverse transformation t^r) appears naturally. Notice that in practical appli- 
cations, both the original Euler energy E and mapped energy -Emap(T) are nonzero. 

We reconstruct the original time t{T) by integrating dt/dr = 1/||lj(-, t(r))||oo and using ([T]): 

tir) = r ,, , ] ,,,, dr' = r . /^JS55 dr' . (8) 

^ ^ J, ||a;(-,t(rO)||oo Jo ^ E ^ ^ 

This allows us to re-word the BKM criterion for singularity of Euler, in terms of the regular 

solution of the mapped equations (|4])-([5]). Let us define the quantity tend: 



Corollary 3. The solution of the 3D Euler system (Tip has a finite-time singularity if and only if 
tend < oo, in which case the singularity is attained at time t = tend in the original variables. D 

This represents a 'duality' between the mapped regular system and the original system: a fast 
decay in time of the mapped field is equivalent to a finite-time blowup of the original field. 



III. PROPOSED NUMERICAL SIMULATIONS OF MAPPED SYSTEM AND THE 
FINITE-TIME BLOWUP HYPOTHESIS OF 3D EULER EQUATIONS 

A. Advantages of the new approach 

The numerical integration of the mapped system (|4])-(|5]) is completely independent of the 
Euler system ([I]), and its main advantages are: 

Adv. 1. Since the profile of vorticity modulus is always bounded by 1, we expect a better 
regularity of all fields involved in the computation of the next timestep data: velocity field and 
damping coefficient /3(t). This will reduce numerical errors. 

Adv. 2. As the profile of vorticity becomes more and more pronounced near its maximum (with 
vorticity modulus equal to 1 at the maximum), the value of vorticity modulus away from the 
peak tends to zero. Accordingly, all numerical noise generated away from the peak should be 
damped to zero. This is due to the damping produced by /3(r), defined in equation ([5]). 
Adv. 3. Due to the boundedness of velocity, the Courant condition will produce a timestep 



which is bounded from below (see Section III C for details on the numerical code) 



Adv. 4. The new approach puts in the same footing the two potential types of solutions of Euler 
equations: the ones that are globally regular and the ones that blow up in a finite time. The 
reason being that for each of the two types, the mapped version is regular. This gives a common 
platform to study the effect of initial conditions on the singular character of the solutions. 
Adv. 5. For solutions of Euler equations that are candidates for finite-time singularity, the 
integration of the mapped regular equations will get further in time with better accuracy than 
the direct integration of Euler equations. Once the new numerical codes to solve the mapped 
equations are validated, our method will serve as a consistency test for the old Euler codes. 



B. Initial conditions 

The guiding principle to set up initial conditions will be to use discrete mirror symmetries, 
which enable us to reach high resolution by saving memory and computation time. Types of 
initial conditions to be studied range from Taylor-Green-like to anti-parallel vortices, and we will 
apply to all these initial conditions a growth-optimisation procedure which is under development. 
First we set up the initial conditions for the field u(x, t) at t = 0: u(x, 0) = uo(x) such that 
V ■ uo(x) = Vx G fi . We construct the initial vorticity: u^oix) = V x uo(x) , and compute 
its supremum norm and the position Y(0) where it attains its maximum value: ||ti;o(-)||oo = 
supxgn |<^o(x)| = |ci;o(Y(0))| . Finally, we compute the initial energy E = ^ J^ |uo(x)p cPx . 



Following the mapping described in Section |ll C[ we now construct the initial conditions 
corresponding to the system (|4])-([5]). Notice that there is no need to explicitly 'solve' for r as 
a function of t. We merely set r(0) = 0. The initial conditions for the mapped velocity are 
thus Umap(x, 0) = uo(x)/| |uJo(-) I loo • The mapped vorticity u>rna,p = V x u^ap satisfies initially 
|^map(x,0)| < 1, Vx G Q, and |u;^ap(Y(0),0)| = 1. 



C. Numerical integration of the mapped equations 

We will integrate numerically the mapped regular system (|4])-(|5]). For this we will modify 
extensively two existing Euler codes, in order to develop a new code where at each time step the 
following extra procedures are implemented: 

• an accurate modelling of the nonlocal damping coefficient /3(r), equation (|5]), using a high- 
order interpolation method to find the instantaneous position Y(r) of the maximum of vorticity 



modulus and then computing the interpolated gradient of velocity at that point; 

• a normalisation algorithm that guarantees that the vorticity modulus be exactly equal to 1 at 
X = Y(r), as it should be analytically, and consequently |ajmap(x, r)| < 1 Vx G i7. 

The former procedure takes a reasonably small fraction of the computation time, depending on 
the required interpolation accuracy. The latter procedure re-scales the fields uniformly in space, 
so it will increase only slightly the computation time. The new code will be pseudospectral and 
parallel, using a Message-Passing-Interface protocol (MP!). Time stepping will follow a high-order 
Runge-Kutta scheme, with a Courant condition to choose the value of the time step. 

The following quantities need to be saved for post-processing analyses that can lead to 
conclusions about the finite-time singularity hypothesis: damping coefficient (|5]), mapped energy 
(6) and mapped circulations: crj(r) = ^^ Umap(r,r) ■ dr , j = 1, . . . ,N , where {Cj}jLi are 
selected fixed contours that depend on the type of mirror symmetry used. 

A reliability time t^ci will be determined for each numerical simulation performed, so that the 
simulation is reliable for times r < Trd. Reliability checks to be produced are: conservation of 
quantities Kj = aj (t) / [Ejaa.p{T)Y^'^ defined in terms of mapped energy and circulation, classical 
resolution studies and analyticity-strip convergence studies [55] . 



D. Searching for late-time trends of singular/nonsingular behaviour 

We proceed to interpolate the saved quantities as functions of time r between and t^cI, i" 
particular the supremum norm of the original vorticity ||aj(-, t(r))||oo = (-E/-Emap ('*")) • 

Notice that in the usual setting of Euler equations, researchers have tried various types of fits 
for the late-time trend of the vorticity supremum norm as a function of time t. For example, a dou- 
ble exponential ||aj(-, t)||oo ~ a exp(exp(6t)) in [51j and a power law ||lj(-, t)||oo ^ c{T^: — tY'^ 
in [55j. We can classify all types of fits in two classes: in one class, the fits that do not represent 
a finite-time singularity and in the other class, the fits that represent explicitly a finite-time 
singularity. In both classes, ||li;(-, t)||oo increases monotonically with time t. 

In contrast, in the new setting in terms of time r and mapped variables, the frontier between 
the two classes is more subtle. According to Theorem 2, the mapped fields are regular for all 
times r < oo. In the two classes, any fit for r — t- oo trends of ||aj(-,t(r))||oo must be continuous. 

In the new setting, there is a more robust means to conclude on the finite-time singu- 
larity hypothesis. From Corollary 3, the blowup time is tend = /o°° ||<-^("j^(''"'))lloo^ (^'^' = 
j^ (-Emap (''"') /-^)^ '^'^' ■ Therefore, a late-r fit for E^^^ir) such that this integral converges, 
provides robust and reliable evidence in favour of the finite-time singularity hypothesis. Con- 
versely, if the integral diverges, the evidence will be against the hypothesis. 

To illustrate these ideas, let us suppose for example that we have obtained a late-time trend 
||ci;(-, t(r))||oo ^ K T^ , where i^ > and 7 > are fit parameters. Using equation rtsl) we solve 
for t(r) and, subsequently, find ||ci;(-, t)||oo by inversion and substitution. We get: 



{l-y)K 
K 



^lnr + to,7 = l llUJIloo <^ KexpKit-to),^ = l 



where to is a constant of integration and i^ > is a constant. We see that the simple power-law 
fit ||a;(-, t(T))||oo ^ K t"' , leads to a prediction of finite-time singularity if and only if the 
exponent 7 > 1 , which is in complete agreement with Corollary 3. In the case 7 > 1, the 
hypothetical singularity time comes directly from equation ([9]), and is not a fit parameter. 



IV. OTHER EQUATIONS THAT ADMIT A MAPPING TO A REGULAR SYSTEM: THE 
CASE OF IDEAL MAGNETOHYDRODYNAMICS 

The analysis in the previous two Sections can be extended to include other equations of 
physical and mathematical interest such as 3D Navier-Stokes, ideal 2D and 3D magnetohydro- 
dynamics (MHD) and even ID inviscid Burgers (see Section [V]). The basic ingredient needed is 
a BKM-type of theorem. For simplicity we discuss briefly the method for the equations of ideal 
2D MHD: 

{dt + u-V)uj = h-Vj, {dt + u-V)4j = 0, (10) 

where V denotes the 2D gradient, the vector fields b (magnetic induction) and u (velocity) are 
defined in terms of the scalars ^ (magnetic potential) and cu (vorticity) by b = V"'"V^ , u = 
V"'"A~^ci;, and the current j is given by j = A^. The perpendicular gradient is defined by 
V"*" = k X V , where k denotes the out-of-plane unit vector. 

A theorem of BKM type, valid for 2D and 3D MHD, was found by Caflisch, Klapper and 
Steele [57j and it states, for the 2D case: 



BKM theorem for MHD. (Caflisch et al.) The solution of the MHD equations (10) is regular 
up to and including time t if and only if the following integral is bounded: 

r{t)= r(||a;(-,t')||oo + ||j(-,Olloo) dt' . U (11) 

JQ 

Corresponding to this time transformation we define mapped fields as follows: 

0(x,r)= „ , ,!^^'''^}., ,„ , ^(x,r)= „ , J^""'^}., ,„ , (12) 

^ ' \\uj{-Moo + \\j{-Moo ^ ^ ||u;(-,t)|U + ||j(-,t)||oo ^ ^ 

along with the mapped velocity U = V"'"A~^f2 and mapped current J = A^ . 

In the new variables we have the identity ||r2(-, r)||oo + ||-^(-7 t)||oo = 1 Vr, which assures 
the boundedness of the supremum norms of mapped vorticity and current. It will be useful to 
define the respective spatial positions Yi(r), Y2(r) where the mapped vorticity il and current 
J attain their respective maxima. We write 

(13) 
(14) 

As discussed in the 3D Euler case for the position of maximum vorticity, the positions Yi(r), Y2(r) 
have discontinuities when competing maxima of vorticity and current coincide. However, in a 
fast-growth scenario, we expect these discontinuities to be distributed as a finite number of iso- 
lated points in the time axis, so after a transient we can consider these functions as continuous. 
The discontinuous case can be treated without difficulty but for simplicity of presentation we 
assume that these functions are continuous. 



m-,r)u = 


= fi(Yi(r),r) , 


Vf](Yi(r),r) = 0, 


^(■,r)U 


= J(j2{r),r) , 


VJ(Y2(r),r) = 0. 



After the bijective mapping (11)-(12), the new variables satisfy the following system of 
equations: 

(9, + U ■ V) fi = B • V J - 7(r) n, (a, + U ■ V) ^ = -^{t) ^ , (15) 

where the damping/amplifying coefficient 7(r) is a quadratic function of the new variables. 
Explicitly, we have 

7(r) = [sgnfiB • VJ]^ + [sgnJB ■ Vn]^ + [2epgdrBpdrUg]^ , (16) 



where sgn represents the signum function, Spg are the components of the 2-dimensional fully 
antisymmetric tensor with £12 = 1, Einstein sum convention over repeated indices is assumed, 
and [F(x,r)]i = F(Yi(r),r), [F(x,r)]2 = F(Y2(r),r), for any function F(x,r). The 
coefficient 7(r) will be damping if and only if /(|U(x, r)p + |B(x, r)p) (Px grows in time. 

In a way completely analogous to the 3D Euler case we can prove the global regularity of the 
MHD mapped system: 



Theorem 4. Assuming smooth initial conditions, the solutions to the mapped system { 15)-(16) 
is regular for all times r > 0. D 

Also, it is possible to reword the BKM-type of theorem for MHD in terms of the mapped 
fields. Defining the original energy (constant) and the mapped energy (not constant) by 

£^^ y"(|u(x,t)p + |b(x,t)n d'x, S^,^{t) ^ ^ y"(|U(x,r)p + |B(x,r)n d'x , (17) 

we obtain ||a;(-, t(r))||oo + ||J(-,^(t))||oo = \/^7^map07. where t(r) denotes the inverse of 
transformation (11). Defining the singularity time tend = /q^ \/£maJp(j)7£ dr , we have: 



Corollary 5. The solution of the MHD system (10) has a finite-time singularity if and only if 



tend < 00, in which case the singularity is attained at time t = tend '" the original variables. D 

Notice that the amplifying coefficient is generically quadratic in the mapped MHD fields, in 
contrast to the Euler case where the coefficient is linear in the fields (for inviscid Burgers, see 
Section |y] the coefficient turns out to be a constant!). This extra nonlinearity might determine 
an extra dealiasing factor in numerical applications, but this is to be established in a forthcoming 
paper. 



V. PROOF OF CONCEPT: AN ILLUSTRATIVE NUMERICAL EXAMPLE 
A. The inviscid Burgers equation in one dimension 

To illustrate the usefulness of the mapping ([S]), we have produced numerical simulations of 
a well-known nonlinear wave equation, called "inviscid Burgers equation" in modern literature, 
although it was introduced and solved at least as early as 1808 by Poisson |58] . 

We have chosen this equation for illustration purposes for three reasons: (i) its solutions are 
known to blow up in a finite time, (ii) it posseses a BKM type of criterion for blowup, and 
(iii) all relevant norms (supremum of gradient, Sobolev) can be found analytically in terms of 
the initial conditions, by the method of characteristics. Point (i) ensures that the example is 
nontrivial because the numerical method is challenging. Point (ii) allows us to apply the time 
mapping in complete analogy to the case of 3D Euler equations. Point (iii) allows us to use 
the analytical expressions for the relevant norms, to produce reliable comparisons and validation 
tests of numerical datasets from both original equation and mapped regular system. 

We consider a spatially periodic setting for a scalar field u{x,t) E M depending on x G M 
and t G [0,tend). such that u{x + 27r,t) = u{x,t), with initial condition u{x,0) = Uo{x) also 
periodic and of class C°°. The field u{x,t) satisfies the differential equation 

ut + uux = 0, (18) 

where the subscript denotes differentiation with respect to the corresponding argument. The 
blowup time tend > is given in terms of the initial condition uo(x) through the well-known 
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formula tend = —i^^xo€[o,2'K]u'Q{xo)y^. Relevant constants of motion in this case are the 
energy E = ^ J^^"" |u(x,t)p dx and 'circulation' ctq = Jq"" u{x,t) dx . 

The analogue of supremum norm of vorticity in 3D Euler equations is for inviscid Burgers the 
supremum norm of m^,, defined by 

\M-,t)\\^ = sup \u^{x,t)\ . (19) 

xe[0,27r] 

For the sake of simplicity of presentation we avoid the situation of competing local maxima 
of the gradients, by assuming that ||^o(-)lloo_~ \''^^^xoe[o,2Tr]UQ{xo)\ . Then, there is a wel 



known analytical expression for the norm (19), obtained by solving equation (18) along the 
characteristics: 

\M;t)\L = r^f (20) 

''cnd f- 

For times t < tend, it can be shown by explicit computations that the solution remains in the class 
C°°, while for t = tend all supremum norms of the spatial derivatives diverge. By direct inspec- 
tion it is therefore evident that, for inviscid Burgers, the analogue of 3D Euler's BKM theorem is: 



Theorem 6. The solution u{x, t) of equation (18) is bounded up to and including time t if and 
only if the following integral is bounded: 



T{t)= / ||«..(-,t')|L dt'.D (21) 

Jo 

The above equation defines the mapped time T{t), and accordingly we define the mapped 

velocity field v{x,t) = u{x,t)/\\ux{-,t)\\oo ■ The bijectively-mapped equation fort>(x, r) is 

Vr + vvx = -a{T) V , (22) 

where the damping/amplifying coefficient is defined by a(r) = —Vx(Y(T),r), and Y(r) is 
the spatial position of the maximum of \vx{x,t)\ . By virtue of the above transformation we 



have |a(r)| = \vx{Y{t),t)\ = \\vx{-,t)\\^ = 1. Therefore, in our mapped equation (22), 
the damping/amplifying coefficient can generally take only two values: either a(r) = 1 or 
a(r) = —1. Discontinuities of this coefficient can arise when two local maxima of \vx{x,t)\ 
compete. To avoid this situation, we have chosen for simplicity the initial condition such that 
IIKlOlloo ^ I i^iiia;oe[o,27r] ^ol^^o)! • So, in our example the damping coefficient turns out to be a 
constant, equal to unity. 

B. Mapping fields 



From equations (20), (21) we obtain 



r(t) = -ln(tend-t), \\ux{-,t{T))\\^ = - — exp T . (23) 

''end 

In numerical simulations of the mapped equation ( [22) ) we have direct access to the mapped 
energy and circulation 

1 /»27r /'27r 

-Eniap(T) = - / |t;(a;,r)|^ da;, a{T) = v{x,t) dx . 

^ Jo Jo 



In a manner completely analogous to what was done in Section |ll D[ we obtain the iden- 
tities ||Mx(-,t(r))||^ = J eJ!^{t) ' ^end = J^ y '^E^''^' dr , and the constant of motion 

J^^ _ ^(t) _ £0 

Correspondingly, the mapped energy and circulations can be obtained explicitly: 

-Emap(i") = (tend)^ E exp(-2r) , a^r) = o-Q^cnd exp(-r) . (24) 
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C. Numerical simulations of inviscid Burgers and mapped system: comparison of data with 

the analytical solution 

For simplicity we will take the initial condition uo{xo) = 1+cosxo , which satisfies the assump- 
tions at the beginning of this Section. The blowup time is tend = 1- The constants of motion for 
inviscid Burgers studied here, energy and circulation, are respectively E = | /q'^ |Mo(a^)P dx = 

l.Svr, o"o = Jq^ uo{x) dx = 271 . We have put deliberately an additive constant in this initial 
condition so that the position of maximum of \ux{x,t)\ moves in time, which is a necessary 
test for the algorithm that searches for the accurate position of the maximum in the numerical 
solution of the mapped equations. 

In order to produce a proper comparison of the data arising from the numerical simulation 



of the inviscid Burgers equation (18) and the mapped equation (22), we approach the problem 
using the following principles: 

• The basic numerical method used for both systems should be identical. We choose as a 
method the classical pseudospectral with leap-frog time advancement (second order accurate), 
which uses a constant time stepping determined by a CFL condition with a typical Courant 
number between 1/5 and 1/20. This method has been used previously in numerical simulations 
of inviscid Burgers [56J. Our method differs slightly with the one used in the cited reference, in 
that we do not mix odd- and even-timestep data every 20 steps. We found that it is not needed 
to do this mixing to obtain convergence. The damping term in the regular system is treated 
using a Crank-Nicholson scheme, which keeps the second-order accuracy. 

• Unlike the integration of the inviscid Burgers equation, the integration of the regular system 
requires a normalisation of the fields every time step, so that the supremum norm is always 
equal to one. For this, at each time step the spatial position of the maximum of \vx{x,t)\ is 
located using an accurate method of order 16, then the function |t>^.(a;,r)| is interpolated at 
that position to obtain ||fa(-, r)||oo, which is used to rescale the field v{x,t). This procedure 
takes only a small fraction (less than 10 percent) of the full computer time used in a time 
stepping. To compensate for this extra time, we have proportionally reduced the CFL number in 
the integration of inviscid Burgers equation when producing comparisons, so that we will always 
be comparing runs with the same spatial resolution, that took approximately the same computer 
time and number of timestep iterations. 

We will compare the numerical data with the analytical solution. We map the data for the 



maximum vorticity obtained from the integration of inviscid Burgers equation (18), using the map 
r = — ln(l —t), so that in all plots below we have r as the time variable. From equation (23) we 
should have ||wa;(-, t(r))||^ = expr. Therefore we plot the error ei(r) = 2(r — In ||^a('!i(T))||oo) 
versus r, where the data used comes from the integration of inviscid Burgers equation. If the 
numerical data is sensible, this error should be close to zero. Figure [T] left panel, shows a 
resolution study with resolutions N = 1024,2048,4096 and 16384. 

For comparison, we use the data for the mapped energy EmapiT), obtained from the inde- 



pendent integration of the mapped system (22). From equation (24) we should have -EmaplT) = 
l.Svr exp(— 2r). We plot, versus r, the error e2{r) = 2r + \n Eraa.p{T) — ln(1.5 7r) which an- 
alytically should be equal to ei(r) and therefore close to zero. Figure [I] right panel, shows a 
resolution study with resolutions A^ = 1024,2048,4096 and 16384. 

Notice that the scale used for both figures is the same. We observe that the convergence 
towards the analytical solution, as one increases the resolution, is much faster for the mapped 
system data than for the original inviscid Burgers system. Also, it is evident that the numerical 
simulation of the mapped system captures the singularity with better accuracy in the time domain 
for a given spatial resolution: figure [2] left panel, shows a comparison of A^ = 16384 data from 
both methods. The error of the numerical data of the mapped system is one order of magnitude 



12 



ei 



10" 



5.x 10 





FIG. 1. Colour online. Left panel: Resolution study of data from integration of inviscid Burgers 
equation. Plots, for different spatial resolutions, of the error ei with respect to the analytical 
solution. Symbols (blue online): '+', 'x', '•' and 'A' correspond respectively to resolutions 1024, 
2048, 4096 and 16384. For reference, the reliability time computed from data of integration 
of mapped system at resolution 16384 is indicated by an arrow (orange online). Right panel: 
Resolution study of data from integration of regular mapped system. Plots, for different spatial 
resolutions, of the error 62 with respect to the analytical solution. Symbols (red online): '+', 'x', 
'o' and 'A' correspond respectively to resolutions 1024, 2048, 4096 and 16384. 

less than the error of the numerical data of the original system. This effect is due in part to the 
fact that the integration of the original inviscid Burgers system uses a constant timestep and 
therefore is bound to misrepresent the solution at times that are close to the singularity time 
4nd = 1- In contrast, the integration of the mapped system also uses a constant timestep but 
the singularity is at Tend = 00 , so the solution is well-resolved as long as the spatial resolution is 
not lost due to the localisation of the structures. 
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FIG. 2. Colour online. Left Panel: Comparison of errors with respect to the analytical solution, 
of datasets obtained from integration of inviscid Burgers equation (error ei, blue bullets) and 
regular mapped system (error 62, red open circles), at fixed resolution 16384. Right panel: Log- 
lin plot of energy spectrum E^ oc \v{k,t)\'^ (orange bullets), of the solution v{x,t) of the regular 
mapped equations at time r = 4.2, using data from numerical integration of regular mapped 
system. The resolution is 16384, corresponding to fcmax = 5461 by the 2/3-rd dealiasing rule. 
It is possible to observe at least two decades of exponential decay of the spectrum (rectangular 
region, light blue online), in agreement with the criterion of reliability commonly used in the 
analiticity-strip method. 



Finally, notice that for a given spatial resolution both original system and regular system have 
a 'reliability time' t^c\ after which the fields are not well-resolved spatially anymore. One way 
to quantify the reliability time is by using the analiticity-strip method [56J, see also [321 |55] 
where the method is applied to 3D Euler equations. At resolution 16384, a calculation using 
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the analiticity-strip method would give Xrei ~ 4.2, corresponding to original time trci ~ -99 (see 
figure 12} right panel, and its caption). An ad-hoc and more conservative way to visualise the 
reliability time is to look at figure[T] right panel. The reliability time for a given spatial resolution 
is approximately the time when the local minimum of the error occurs. For example, at resolution 
16384 the reliability time is Trd ~ 3.78, corresponding to original time tj-d ~ 0.98 . It is clear from 
figure |2} left panel, that the data from the mapped system agrees with the analytical solution 
up to times very close to the corresponding reliability time, while the data for inviscid Burgers 
system begins to diverge earlier. 

As an extra check, using data from the numerical solution of the mapped regular system, we 
have confirmed the constancy of K = (t{t)/ (-EmaplT)) = ctq/ E^^"^ = 2 (7r/3) ' for all runs, 
within one part in 10^^ for resolution 1024, and well beyond the reliability time. 



VI. CONCLUSION AND DISCUSSION 

We have established, for several equations of physical interest, a nonlinear time transformation 
and field rescaling that is based on a BKM-type of theorem, and which maps bijectively the 
original equations to a system that is globally regular in time. 

We have learnt from the simple example of inviscid Burgers that the numerical integration 
of a mapped regular system produces more accurate and reliable results than the integration 
of the original system. The only way that the integration of the original system can match the 
regular system, is by upgrading the numerical method to an adaptive time stepping, based on the 



transformation (21), so that dr = dt x \\ux{-,t)\\ao is kept constant. In addition, we validated 
a procedure of spatial interpolation (of order 16) used to compute accurately the value of the 
supremum norm ||fx(-, t)||oo- This procedure is essential for the numerical integration of the 
corresponding mapped regular system, not only in inviscid Burgers but also in 3D Euler, Navier 
Stokes and magnetohydrodynamics. It is worth mentioning that we also used this interpolation 
procedure to compute ||Mx(-,i)||oo from the numerical data of inviscid Burgers equations. If we 
had just computed the collocation-point maximum value of the same quantity, we would have 
obtained a prohibitive amount of spurious oscillations (figure not shown). 

In a forthcoming paper the method will be implemented to integrate the 3D Euler equations. 
We expect the integration of the regular mapped system to produce a significant improvement 
in accuracy. However in this case no analytical solution is at hand, so in order to draw robust 
conclusions on the trends we will need to make a late-time fit along with a careful determination 
of the reliability time. We are considering the possibility to include also a time-dependent spatial 
deformation in order to capture the localised structures as they develop. This procedure is in 
the spirit of l59i ;60J and the recent review [61j . 
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